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Abstract. In this paper, we investigate the efficiency of various strategies for subdividing 
polynomial triangular surface patches. We give a simple algorithm performing a regular 
subdivision in four calls to the standard de Casteljau algorithm (in its subdivision version). 
A naive version uses twelve calls. We also show that any method for obtaining a regular 
subdivision using the standard de Casteljau algorithm requires at least 4 calls. Thus, our 
method is optimal. We give another subdivision algorithm using only three calls to the 
de Casteljau algorithm. Instead of being regular, the subdivision pattern is diamond-like. 
Finally, we present a "spider-like" subdivision scheme producing six subtriangles in four calls 
to the de Casteljau algorithm. 
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1 Introduction 



In this paper, we investigate the efficiency of various strategies for subdividing polynomial 
triangular surface patches. Subdivision methods based on a version of the de Casteljau 
algorithm splitting a control net into control subnets (see Farin ^) were investigated by 
Goldman Boehm and Farin j2j, Bohm [4 , and Seidel (see also Boehm, Farin, and 
Kahman and Filip jH]). However, except for Bohm jl], these papers are not particularly 
concerned with minimizing the number of calls to the standard de Casteljau algorithm. 
Furthermore, some of these papers (notably Goldman ^T]) use a version of the de Casteljau 
algorithm computing a 5-dimensional simplex of polar values, which is more expensive than 
the standard 3-dimensional version. In this paper, we give a simple algorithm performing a 
regular subdivision in four calls to the standard de Casteljau algorithm (in its subdivision 
version). A naive version uses twelve calls. We also show that any method for obtaining a 
regular subdivision using the standard de Casteljau algorithm requires at least 4 calls. Thus, 
our method is optimal. We give another subdivision algorithm using only three calls to the 
de Casteljau algorithm. Instead of being regular, the subdivision pattern is diamond-like. 
Finally, we present a "spider-like" subdivision scheme producing six subtriangles in four 
calls to the de Casteljau algorithm. Some familiarity with affine spaces and affine maps is 
assumed. Details can be found in Farin [3, Berger Ppj, or Gallier 0. 

2 The Polar Form Approach to Polynomial Triangular 
Surface Patches 

The deep reason why polynomial triangular surface patches can be effectively handled in 
terms of control points is that multivariate polynomials arise from multiaffine symmetric 
maps (see Ramshaw [13], Farin [3 IH], Hoschek and Lasser ^21, or Gallier |^). Denoting 
the affine plane as V, traditionally, a polynomial surface in R" is a function F:V M", 
defined such that 

Xi = Fi(M,f), 
Xn = Fn{u,v), 

for all u,v E M., where Fi{U,V), . . . , Fm{U, V) are polynomials in M[f/, V]. Given a natural 
number m, if each polynomial Fi{U, V) has total degree < m, we say that F is a polynomial 
surface of total degree m. The trace of the surface F is the set F{V). 

Now, given a polynomial surface F of total degree m in some affine space S (typically 
M'^), there is unique symmetric and multiaffine map /: V"^ — > S such that 

F{u,v) = f{{u,v),...,{u,v)), 
^ ' 
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for all {u, v) G V. The symmetric and multiaffine map / associated with F is called the 
polar form of F. 

The above result is not hard to prove. Using linearity, it is enough to deal with a single 
monomial. Given a monomial U^V^, with h + k = d < m, it is easily shown that the 
symmetric multiaffine / form corresponding to U^V^ is given 

IUJC{l,...,m} i&I jeJ 

inJ=0 

card(/)=/i, card( J)=fc 

Recall that a map /: M'^ — > is affine if 

/((l_A)a + A6) = (l-A)/(a) + A/(6), 
for all a, be M"', and all A G M. A map /: R*^ x ■ ■ ■ x R*^^ -> R" is multiaffine if it is affine in 

m 

each of its arguments. A map /: R*^ x ■ — x R'^^ — > R" is symmetric if it does not depend on 

m 

the order of its arguments, i.e., /(a^(i), . . . , a-K{m)) = /(o-i, • • • , dm), for all ai, . . . , am, and all 
permutations vr. 

As an example, consider the following surface known as Enneper's surface: 

F^{U,V) = U-^ + UV^ 

F,{U,V) = V-^ + U''V 
F^{U,V) = f/2 - 1/2. 

We get the polar forms 

/l((f/l,ri),(f/2,V2),([/3,l^3)) 



U1 + U2 + U3 U1U2U, 
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+ 



3 3 
U1V2V3 + U2ViV3 + U-iViV2 



f2m,V,),{U2,V2),{Us,V,)) 



f3m,V,),iU2,V2),iUs,Vs)) 



3 

Vi + V2 + V3 V1V2V3 



+ 



3 3 
UXU2V3 + U1U3V2 + U2U3V1 



3 

U1U2 + U1U3 + U2U3 V^V2 + + V2V3 



Furthermore, it turns out that any symmetric multiaffine map /: — £^ is uniquely 
determined by a family of ('"+iK"^+^) points (where 8 is any affine space, say R"). Let 

A™ = {(z, J, fc) G I z + J + = m}. 

The following lemma is easily shown (see Ramshaw jT3] or Gallier P). 
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Lemma 2.1 Given an affine frame Arst in the plane V, given a family {bij^k)(i,j,k)eAm c"/ 
(m+iKm+2) p^j^^^g there is a unique surface F:V ^ £ of total degree m, defined by a 

symmetric m-affine polar form f: — > S, such that 

i j k 

for all {i,j,k) e A^. Furthermore, f is given by the expression 

^kjf i ^i ' ^' ' g, ■ ^. , S , t, . ,t^ , 

I,J,K disjoint 

where ai — Xir + /liS + Uit, with \i + + Ui — 1, and 1 <i <m. 

For example, with respect to the standard frame Arst = ((1, 0, 0), (0, 1, 0), (0, 0, 1)), we 
obtain the following 10 control points for the Enneper surface: 



/(r, r, r) 
(3. 0.1) 

f{r,r,t) f{r,r,s) 

2 1 2 2 1 

^3'°'3^ ^3'3'3^ 

f{r,t,t) f{r,s,t) f{r,s,s) 

(i.0,0) (iio) 

f(t,t,t) J(s,t,t) /(■':«■«) f{s,s,s) 

("•"•o) (o.i.o) (o,|-i) (o,|-i) 

A family Af = {bi,j^k)(i,j,k)eAm of ("^+^)^(™+^) points in S is called a (triangular) control 
net, or Bezier net. Note that the points in 

= {(i, j,k)eN^ \ i+j + k^ m}, 

can be thought of as a triangular grid of points in V. For example, when m = 5, we have 
the following grid of 21 points: 









500 












401 




410 








302 




311 




320 




203 




212 




221 




230 


104 


113 




122 




131 


140 


005 014 




023 




032 




041 050 
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We intentionally let i be the row index, starting from the left lower corner, and j be the 
column index, also starting from the left lower corner. The control net Af = {bij^k)(i,j,k)eAm 
can be viewed as an image of the triangular grid in the affine space S. It follows from 
Lemma ITT] that there is a bijection between polynomial surfaces of degree m and control nets 
(jj,fc)eAm- It should also be noted that there are efficient methods for computing 
control nets from parametric definitions, but this will be published elsewhere. 

In the next section, we review a beautiful algorithm to compute a point F{a) on a surface 
patch using affine interpolation steps, the de Casteljau algorithm. 



3 The de Casteljau Algorithm for Triangular Patches 

In this section, we explain in detail how the de Casteljau algorithm can be used to subdivide a 
triangular patch into three subpatches. For more details, see Farin [711^1, Hoschek and Lasser 
jl2j . Risler ^3], or Gallier [H]. In the next section, we will use versions of this algorithm to 
obtain a triangulation of a surface patch using recursive subdivision. 

Given an affine frame Arst, given a triangular control net J\f = {bij^k){i,j,k)£A^, recall 
that in terms of the polar form /: V"^ £ of the polynomial surface F:V ^ £ defined by 
jV, for every (i, j, k) G A^,, we have 

bi,j,k = / ( t", . , r , S, . .J , 5 , t, . ^. , t ) . 

i j k 

Given a = Xr + fis + ut in V, where X + fi + u = 1, in order to compute F{a) = f{a, . . . , a), the 
computation builds a sort of tetrahedron consisting of m + 1 layers. The base layer consists 
of the original control points in Af, which are also denoted as j_fc)(i,i,fc)eA™- The other 
layers are computed in m stages, where at stage I, 1 < I < m, the points j,fc)(j,j,fc)gA„_, 
are computed such that 

^i,j,k = ^b\^l,j,k + /^^IZ+Lfc + ^^\,lk+l- 

During the last stage, the single point b^Q^ q is computed. An easy induction shows that 

bij,k = f r^_^^, s^_^^^, t^^^^) , 

I i j k 

where k) G Am-i, and thus, F{a) = b'^o^Q- 

Assuming that a is not on one of the edges of Arst, the crux of the subdivision method 
is that the three other faces of the tetrahedron of polar values b\ j f^ besides the face corre- 
sponding to the original control net, yield three control nets 



Mast = (6o j. 



j,k){i,j,k)eeAm^ 
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corresponding to the base triangle Aast, 

Nrat = (&-,o,ik)w,fc)eeA^, 
corresponding to the base triangle Arat, and 

corresponding to the base triangle Arsa. If a belongs to one of the edges, say rs, then the 
triangle Arsa is flat, i.e. Arsa is not an aflne frame, and the net Afrsa does not define the 
surface, but instead a curve. However, in such cases, the degenerate net Afrsa is not needed 
anyway. 

From an implementation point of view, we found it convenient to assume that a triangular 
net J\f — {bi,j,k){i,j,k)eAm represented as the list consisting of the concatenation of the m + 1 
rows 

i.e., 

/(tV^^, /(^V^^, s, t^^_^^, . . . , f(r\^^^, s^_^^, t), f(r, . . . ,r, s, . . . , s), 

i m—i i m—i—1 i m—i—1 

where < i < m. As a triangle, the net M is hsted (from top-down) as 



m—l m—1 m 



m—l m—l 

fir, . . . ,r) 



The main advantage of this representation is that we can view the net A/" as a two- 
dimensional array net, such that net[i,j] — bi,j,k (with i -\- j -\- k — m). In fact, only a 
triangular portion of this array is filled. This way of representing control nets fits well with 
the convention that the affine frame Arsi is represented as follows: 
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s 



Figure 1: An affine frame 

Instead of simply computing F{a) = 6™q q, the de Casteljau algorithm can be easily 
adapted to output the three nets Mast, Mrat, and Mrsa. The function sdecasS does that. 
We also found it convenient to write three distinct functions subdecasSra, subdecasSsa, 
and subdecasSta, computing the control nets with respect to the affine frames Aast, Aart, 
and Aars. An implementation in Mathematica can be found in Gallier [Hj. 

4 Regular Subdivision Of Triangular Patches 

If we want to render a triangular surface patch F defined over the affine frame Arst, it seems 
natural to subdivide Arst into the three subtriangles Aars, Aast, and Aart, where a = 
(1/3, 1/3, 1/3) is the center of gravity of the triangle Arst, getting new control nets Afars, 
Mast and Mart using the functions described earlier, and repeat this process recursively. 
However, this process does not yield a good triangulation of the surface patch, because no 
progress is made on the edges rs, st, and tr, and thus, such a triangulation does not converge 
to the surface patch. Thus, in order to compute triangulations that converge to the surface 
patch, we need to subdivide the triangle Arst in such a way that the edges of the affine 
frame are subdivided. There are many ways of performing such subdivisions, and we will 
propose a method which has the advantage of yielding a very regular triangulation, and of 
being very efficient. In fact, we give an optimal method for subdividing an affine frame using 
four calls to the standard de Casteljau algorithm in its subdivision version. A naive method 
would require twelve calls. 

Goldman ^T] proposed several subdivision algorithms, including one for splitting a trian- 
gular patch into four triangular subpatches, but his methods use a generahzed version of the 
de Casteljau algorithm computing a 5-simplex of polar values. These methods are illustrated 
graphically in Boehm and Farin |2]. It should be noted that Boehm and Farin do mention 
that it is possible to compute the control net w.r.t. a new affine frame from the control 
net w.r.t. an original affine frame in three calls to the standard de Casteljau algorithm. 
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However, they do not explain how to spht a triangular patch into four subpatches using four 
calls to the standard de Casteljau algorithm. 

Goldman's subdivision methods can be justified in a very simple way as shown by Sei- 
del pni- Given a surface F of total degree m defined by a triangular control net M = 
{bij,k){i,j,k)eAm^ w.r.t. the affine frame Arst, for any n points Pi = UiV + ViS + Wit (where 
Ui+Vi+Wi = 1), the following (n+2)-simplex of points b^lf^'^"' where i+j + k+li + . . . + = m 
is defined inductively as follows: 

7,0,. ..,0 _ , 
"i,j,k ~ "hj,ky 

~ "'i "i+l,j,k '^^i "i,j+l,k ^fi "i,j,k+i ' 

where 1 < h < n. 

If / is the polar form of F, it is easily shown that 




i j k h In 



For n = 0, F{p) = &ooo; standard de Casteljau algorithm. For n = 3, 

(^o,o!o')(;i,«2,i2)6A™ is a control net of F w.r.t. ApipsPa- 

In particular, if Pi,P2,P3 are chosen on the edges of Arst, the subnets for the four sub- 
patches are obtained. Goldman observes that some of the nets involved in the computation 
are trivial, but still, a 5-simplex of polar values is computed. 

It was brought to our attention by Gerald Farin (and it is mentioned in Remark 2 of 
Seidel's paper [16j, page 580) that Helmut Prautzsch showed in his dissertation (in German) 
[T^ that regular subdivision into four subtriangles can be achieved in four calls to the 
standard de Casteljau algorithm. Prautzsch's method is briefly described in Bohm ^[j, page 
348 (figure) and page 349 (in fact, with a typo, one of the barycentric coordinates listed as 
(|, |, l) should be (|, |,0)). The order in which the four patches are obtained is slightly 
different from ours. Since Prautzsch's algorithm has not been discussed more extensively in 
the literature, we feel justified in presenting our method. 

The subdivision strategy that we will follow is to divide the affine frame Arst into four 
subtriangles Aabt, Abac, Acrb, and Asca, where a = (0,1/2,1/2), b = (1/2,0,1/2), and 
c = (1/2, 1/2,0), are the middle points of the sides st, rt and rs respectively, as shown in 
the diagram below: 
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Figure 2: Subdividing an affine frame Arst 



First, wc show that any method using the standard version of the de Casteljau algorithm 
for subdividing a triangular patch into 4 subpatches forming a regular pattern as above 
requires 4 calls. The crux of the argument is that a call to the de Casteljau algorithm in 
its subdivision version produces three subpatches containing only one new corner. We want 
to produce the four subpatches abt, crb, sea, and abc. After one subdivision step, we have 
three patches each involving exactly one of a, b, c. After two subdivision steps, we have six 
subspatches only two of which involve exactly two of a, b, c, since we can only subdivide 
a single patch, and since this patch only has one of a, b, c. Thus, at least three steps are 
needed to produce four subpatches involving at least two of a, b, c. If we produced the patch 
abc during the third subdivision step, we would have three patches involving exactly two 
of a, b, c, but the subdivision step that produced abc also produces two patches sharing the 
same vertex from (r, s,t). However, abt, crb, and sea do not share a vertex from {r,s,t). If 
abc was not produced during the third step, at least four steps are needed. Therefore, in all 
cases, at least four steps are needed to produce the required four subpatches. 

We now present our algorithm. The first step is to compute the control net for the afiine 
frame Abat. This can be done using two steps. In the first step, split the triangle Arst into 
the two triangles Aart and Aars, where a — (0, 1/2, 1/2) is the middle of st. Using the 
function sdecas3 (with a = (0,1/2,1/2)), the nets Mart, Afast, and Mars arc obtained, 
and we throw away Mast (which is degenerate anyway). Then, we split Aart into the two 
triangles Abat and Abar. For this, we need the barycentric coordinates of b with respect to 
the triangle Aart, which turns out (0, 1/2, 1/2). Using the function sdecasS, the nets Mbat, 
Mbrt, and Mbar are obtained, and we throw away Mbrt. 
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Figure 3: Computing the nets J\fbat, Nbar and Mars from Mrst 

We will now compute the net Mcas from the net Mars. For this, wc need the barycentric 
coordinates of c with respect to the triangle Aars, which turns out to be (0, 1/2, 1/2). Using 
the function subdecasSsa, the net Mcas is obtained. 




Figure 4: Computing the net Mcas from Mars 



We can now compute the nets Mcbr and Mcba from the net Mbar. For this, we need 
the barycentric coordinates of c with respect to the affine frame Abar which turns out to be 
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(— 1, 1, 1). Using the function sdecasS, the snet Mchr, Near, and Mcha are obtained, and 
we throw away Hear. 




Figure 5: Computing the nets Nchr and Ncha from J\fbar 



Finally, we apply transposej to the net Afbat to get the net Afabt. transposek to Hcba 
to get the net Afbac, transposej followed by transposek to the net Mcbr to get the net 
J\fcrb, and transposek twice to Ncas to get the net J\fsca, 




Figure 6: Subdividing Arst into Aabt, Abac, Acrb, and Asca 
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Thus, using four calls to the de Casteljau algorithm, we obtained the nets Afabt, jVbac, 
Afcrb, and Msca. 

Remarks: 

(1) For debugging purposes, we assigned different colors to the patches corresponding to 
N'abt, Afbac, Afcrb, and Afsca, and we found that they formed a particularly nice 
pattern under this ordering of the vertices of the triangles. In fact, Afabt is blue, Afbac 
is red, Afcrb is green, and Afsca is yellow. 

(2) In the last step of our algorithm, the subdivision step is performed with respect to a 
point of barycentric coordinates (—1, 1, 1). One might worry that such a step involving 
a nonconvex combination is a source of numerical instability. We tested our algorithm 
on many different examples, and so far, without running into any problem. We also 
believe that such a nonconvex step is unavoidable if the standard de Casteljau algorithm 
(building a simplex of polar values of dimension 3) is used, but we are unable to prove 
this. 

The subdivision algorithm just presented has been implemented in Mathematica, see 
Gallier 0. The subdivision method is illustrated by the following example of a cubic patch 
specified by the control net 

net = {{0, 0, 0}, {2, 0, 2}, {4, 0, 2}, {6, 0, 0} , 
il, 2, 2}, {3, 2, 5}, {5, 2, 2}, 
{2, 4, 2}, iA, A, 2}, {3, 6, 0}}; 
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6 

Figure 7: Subdivision, 1 iteration 
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Figure 9: Subdivision, 3 iterations 



After only three subdivision steps, the triangulation approximates the surface patch very 
well. 

As another example of the use of the above functions, we can display a portion of a well 
known surface known as the "monkey saddle" , defined by the equations 

X = u, y = V, z = — Suv'^. 

Note that z is the real part of the complex number [u + iv)^. It is easily shown that the 
monkey saddle is specified by the following triangular control net monknet over the standard 
affine frame Arst, where r = (1, 0, 0), s = (0, 1, 0), and t = (0, 0, 1). 

monknet = {{0, 0, 0}, {0, 1/3, 0}, {0, 2/3, 0}, {0, 1, 0}, 
{1/3, 0, 0}, {1/3, 1/3, 0}, {1/3, 2/3, -1}, 
{2/3, 0, 0}, {2/3, 1/3, 0}, {1, 0, 1}}; 

We actually display the patch over the rectangle [—1, 1] x [—1, 1]. This can be done by 
splitting the square into two triangles, and computing control nets with respect to these 
triangles. This is easy to do, and it is explained for example in Gallier [H]. Subdividing both 
nets 3 times, we get the following picture. 



14 



Figure 10: A monkey saddle, triangular subdivision 

5 A Diamond-Shape Strategy For Subdivision 

The strategy of the previous section was to split the affine frame Arts into four congruent 
subtriangles. We were able to do this using four calls to the de Casteljau algorithm and we 
showed that it is not possible to do it in fewer calls. 

However, it is possible to split the affine frame into four subtriangles using only three 
calls to the de Casteljau algorithm. The method consists in splitting the triangle Arst into 
the four subtriangles Abat, Abar, Acas, and Acar: 
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Figure 11: Diamond-style subdivision of an affine frame Arst 



This can be done by first computing the nets Mart and Mars, which can be done in 
one call to sdecasS (dropping Mast). Next, we split Aart into the two triangles Ahat and 
Abar. For this, we need the barycentric coordinates of b with respect to the triangle Aart, 
which turns out (0, 1/2, 1/2). Using the function sdecasS, the nets Mbat, Mbrt, and Mbar 
are obtained, and we throw away Mbrt. Finally, we split Aars into the two triangles Acas 
and Acar. For this, we need the barycentric coordinates of c with respect to the triangle 
Aars, which turns out (0, 1/2, 1/2). Using the function sdecasS, the nets Mcas, Mcrs, and 
Mcar are obtained, and we throw away Mcrs. 

An implementation of the method is given in Gallier [Oj. 

The result of subdividing two of three times reveals some diamond-shape subdividion 
patterns. For example, after three iterations, the dome surface is subdivided as follows: 
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Figure 12: Diamond-style subdivision, 3 iterations 



6 A Spider- Web Strategy For Subdivision 

Is is also possible to split the affine frame into six subtriangles using only four calls to the 
de Casteljau algorithm. 

s 




Figure 13: spider-web subdivision of an affine frame Arst 
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The triangle Arst is subdivided into the triangles Abgt, Abgr, Aagt, Aags, Acgs, and 
Acgr, as shown above, where g is the center of gravity. This can be done by first computing 
the nets Afgrt, Mgst, and Mgrs, which can be done in one call to sdecasS. We split Agrt 
into the two triangles Ahgt and Ahgr using sdecasS (throwing away Abrt). We split Agst 
into the two triangles Aagt and Aags using sdecasS (throwing away Aast). Finally, we 
split Agrs into the two triangles Acgs and Acgs using sdecasS (throwing away Acrs). 

The result of subdividing recursively yields spider-web like patterns. For example, after 
three iterations, the dome surface is subdivided as follows: 



6 
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Figure 14: Spider-web style subdivision, 3 iterations 

7 Conclusion 

We have presented various strategies for subdividing polynomial triangular surface patches. 
We gave an algorithm performing a regular subdivision in four calls to the standard de 
Casteljau algorithm, and we showed that this method for obtaining a regular subdivision is 
optimal. We gave another subdivision algorithm using only three calls to the de Casteljau 
algorithm. Instead of being regular, the subdivision pattern is diamond-like. Finally, we 
presented a "spider-web" pattern subdivision scheme producing six subtriangles in four calls 
to the de Casteljau algorithm. These methods immediately apply to rational surface patches 
(Gallier [10]). 
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An amusing effect is obtained from tlie regular subdivision scheme if we omit the central 
triangle Abac. We obtain a "fractalized" representation of the surface patch, in the sense 
that a Sierpinski gasket pattern is laid onto the patch! It would be interesting to investigate 
other subdivision strategies and the patterns that they induce, especially if the triangles are 
colored in various recursive manners. 
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